function inertiaRatio = calcInertiaRatio( vecPhi, vecR )
vecX = vecR.*cos(vecPhi);
vecY = vecR.*sin(vecPhi);

N = length(vecPhi);
%% area, center-of-mass
accArea = 0;
accIx = 0;
accIy = 0;

for k = 0:N-1
    l = mod(k+1,N);
    xi = vecX(k+1); xj = vecX(l+1);
    yi = vecY(k+1); yj = vecY(l+1);
    
    accArea = accArea + xi*yj - xj*yi;
    accIx = accIx + (xi*yj - xj*yi)*(xi + xj);
    accIy = accIy + (xi*yj - xj*yi)*(yi + yj);
end
area = 1/2 * accArea;
Ix = 1/6 * accIx;
Iy = 1/6 * accIy;

xb = Ix/area;
yb = Iy/area;
%% moment-of-area tensor
accT11 = 0;
accT22 = 0;
accT12 = 0;
for k = 0:N-1
    l = mod(k+1,N);
    xi = vecX(k+1)-xb; xj = vecX(l+1)-xb;
    yi = vecY(k+1)-yb; yj = vecY(l+1)-yb;
    
    accT11 = accT11 + (xi*yj - xj*yi)*(yi^2 + yi*yj + yj^2);
    accT22 = accT22 + (xi*yj - xj*yi)*(xi^2 + xi*xj + xj^2);
    accT12 = accT12 + (xi*yj - xj*yi)*(2*xi*yi + xi*yj + xj*yi + 2*xj*yj);
end
T11 = 1/12 * accT11;
T22 = 1/12 * accT22;
T12 = 1/24 * accT12;
T = [T11 T12; T12 T22];

[v,d] = eig(T); eigenvals = diag(d); mask = abs(v(1,:)) > abs(v(2,:));
if sum(mask) == 0
    mask = [true false];
end
e(1) = eigenvals(mask); e(2) = eigenvals(not(mask));
inertiaRatio = e(1) / e(2);
end